Feature analysis

1 Setup

library(magrittr)
library(tidyverse)
library(glue)
library(arrow)
library(matric)

batch <- params$batch
futile.logger::flog.info(glue("Batch = {batch}"))

INFO [2021-04-17 12:15:27] Batch = 2020_11_04_CPJUMP1

experiment <- as.data.frame(params$experiment)
futile.logger::flog.info(glue("Experiment = {experiment}"))

INFO [2021-04-17 12:15:27] Experiment = Standard
INFO [2021-04-17 12:15:27] Experiment = Compound
INFO [2021-04-17 12:15:27] Experiment = U2OS
INFO [2021-04-17 12:15:27] Experiment = 48

data_level <- params$data_level
futile.logger::flog.info(glue("Data level = {data_level}"))

INFO [2021-04-17 12:15:27] Data level = normalized

if (params$normalization == "whole_plate") {
  norm_suffix <- ""
} else if (params$normalization == "negcon") {
  norm_suffix <- "_negcon"
}
sprintf("Normalization = %s. Using suffix = '%s'", params$normalization, norm_suffix)

[1] "Normalization = negcon. Using suffix = '_negcon'"

filename_prefix_profiles <- glue("{batch}_all_{data_level}{norm_suffix}")

2 Load

parquet_file_recoded <-
  glue("../collated/{batch}/{filename_prefix_profiles}_recoded_aggregated.parquet")

if (file.exists(parquet_file_recoded)) {
  parquet_file <- parquet_file_recoded
  
  futile.logger::flog.info(glue("Loading {parquet_file} ..."))
  
  stopifnot(file.exists(parquet_file))
  
  profiles <-
    read_parquet(parquet_file)
  
  variables <-
    str_subset(names(profiles), "Metadata_", negate = TRUE)
  
} else {
  
  parquet_file <-
    glue("../collated/{batch}/{filename_prefix_profiles}.parquet")
  
  futile.logger::flog.info(glue("Loading {parquet_file} ..."))
  
  stopifnot(file.exists(parquet_file))
  
  profiles <-
    read_parquet(parquet_file)
  
  if(!is.null(experiment)) {
    profiles <- profiles %>% inner_join(experiment)
  }
  
  variables <-
    str_subset(names(profiles), "Metadata_", negate = TRUE)
  
  metadata_cols <-
    str_subset(names(profiles), "Metadata_")
  
  variables <-
    params$variable_groups %>%
    unlist() %>%
    map(function(pattern)
      str_subset(variables, pattern = pattern)) %>%
    unlist()
  
  profiles <-
    profiles %>%
    select(all_of(c(metadata_cols, variables)))
  
  variables <-
    str_subset(names(profiles), "Metadata_", negate = TRUE)
  
  profiles <-
    profiles %>%
    group_by(across(params$strata_replicate)) %>%
    summarize(across(all_of(variables), median),
              .groups = "keep")
  
  attr(profiles, "variable_groups") <- params$variable_groups
  
  parquet_file <-
    glue("../collated/{batch}/{filename_prefix_profiles}_recoded_aggregated.parquet")
  
  futile.logger::flog.info(glue("Writing {parquet_file_recoded} ..."))
  
  profiles %>%
    write_parquet(parquet_file,
                  compression = "gzip",
                  compression_level = 9)
}

INFO [2021-04-17 12:15:28] Loading ../collated/2020_11_04_CPJUMP1/2020_11_04_CPJUMP1_all_normalized_negcon.parquet ...

Joining, by = c("Metadata_cell_line", "Metadata_timepoint", "Metadata_experiment_type", "Metadata_experiment_condition")

INFO [2021-04-17 12:16:21] Writing ../collated/2020_11_04_CPJUMP1/2020_11_04_CPJUMP1_all_normalized_negcon_recoded_aggregated.parquet ...

3 Transform

3.1 Functions

get_beta <- function(x = 5, y = 0.99) {
  -log(1 / y - 1) / x
}

abs_logistic <- function(x, beta = get_beta()) {
  2 / (1 + exp(-abs(x) * beta)) - 1
}

cytominer_variable_group_enrichment <-
  function(population,
           variables,
           variable_groups,
           sigmoid_function = abs_logistic,
           ...) {
    
    variables_group_lists <-
      variable_groups %>%
      unlist() %>%
      names() %>%
      set_names() %>%
      purrr::map(function(variable_group) {
        stringr::str_subset(variables, variable_group)
      }) 

    population_data_transformed <-
      variables_group_lists %>%
      map_dfc(function(variables_group_list) {
        population %>%
          dplyr::select(all_of(variables_group_list)) %>%
          dplyr::mutate(across(all_of(variables_group_list),
                               sigmoid_function,
                               ...)) %>%
          dplyr::ungroup() %>%
          dplyr::rowwise() %>%
          dplyr::transmute(n_above = round(sum(
            dplyr::c_across(everything()) /
              length(variables_group_list),
            na.rm = T
          ), 2)) %>%
          dplyr::pull(n_above)
      })
    
    population_metadata <-
      population %>%
      dplyr::select(-all_of(variables))
    
    enriched <-
      dplyr::bind_cols(population_metadata,
                       population_data_transformed)
    
    enriched
    
  }

3.2 Execute

profiles_enriched <-
  profiles %>%
  group_by(across(params$strata_replicate)) %>%
  summarize(across(all_of(variables), median), .groups = "keep") %>%
  group_by(across(params$nesting_level_0)) %>%
  summarise(
    cytominer_variable_group_enrichment(
      cur_data_all(),
      variables = variables,
      variable_groups = params$variable_groups
    ),
    .groups = "keep"
  ) %>%
  ungroup()

variables_enriched <-
  str_subset(names(profiles_enriched), "Metadata_", negate = TRUE)

4 Plot

profiles_enriched %>% 
  pivot_longer(-matches("Metadata"), 
               names_to = "feature_group")%>% 
  ggplot(aes(feature_group, value)) + 
  geom_boxplot() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

profiles_enriched %>% 
  pivot_longer(-matches("Metadata"), 
               names_to = "feature_group") %>% 
  ggplot(aes(feature_group, value, group = Metadata_Well)) + 
  geom_line(alpha = .09) + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
  facet_wrap(~Metadata_negcon_or_other, ncol = 1)

library(plotly)

df <- profiles_enriched

dimensions <- list(
  list(label = "Intensity_RNA",
       values = ~ Intensity_RNA),
  list(label = 'Intensity_DNA',
       values = ~ Intensity_DNA)
)

  
df <- df %>% select(all_of(variables_enriched))

dimensions <-
  variables_enriched %>%
  map(function(variable) {
    list(label = variable,
         values = df[[variable]])
  })

df %>%
  plot_ly(width = 1000, 
          height = 600,
          hoverinfo = "text") %>%
  add_trace(type = 'parcoords',
            dimensions = dimensions)

'parcoords' objects don't have these attributes: 'hoverinfo'
Valid attributes include:
'type', 'visible', 'name', 'uid', 'ids', 'customdata', 'meta', 'stream', 'transforms', 'uirevision', 'domain', 'labelangle', 'labelside', 'labelfont', 'tickfont', 'rangefont', 'dimensions', 'line', 'idssrc', 'customdatasrc', 'metasrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'
'parcoords' objects don't have these attributes: 'hoverinfo'
Valid attributes include:
'type', 'visible', 'name', 'uid', 'ids', 'customdata', 'meta', 'stream', 'transforms', 'uirevision', 'domain', 'labelangle', 'labelside', 'labelfont', 'tickfont', 'rangefont', 'dimensions', 'line', 'idssrc', 'customdatasrc', 'metasrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

5 Save

attr(profiles, "variable_groups") <- params$variable_groups

parquet_file <-
  glue("../collated/{batch}/{filename_prefix_profiles}_enriched.parquet")

futile.logger::flog.info(glue("Writing {parquet_file_recoded} ..."))

INFO [2021-04-17 12:17:04] Writing ../collated/2020_11_04_CPJUMP1/2020_11_04_CPJUMP1_all_normalized_negcon_recoded_aggregated.parquet ...

profiles_enriched %>%
  write_parquet(parquet_file,
                compression = "gzip",
                compression_level = 9)